clc;
clear;

basis_function_plot("Non-com", "x");

%% basis_function_plot
function basis_function_plot(basis_type, derivative_type, flag)
[x, y] = meshgrid(0:0.0125:1, 0:0.0125:1);
z = basis_function_ref(x(:), y(:), basis_type, derivative_type);
if 2 == nargin
    Nb = size(z,2);
else
    Nb = min(size(z,2), flag);
end
triangle_region = @(x,y) (x >= 0) & (x <= 1) & (y >= 0) & (y <= 1) & (y <= -x+1);
for j = 1:Nb
    Z = reshape(z(:,j), size(x,1), size(x,2));
    Z(~triangle_region(x, y)) = NaN;
    surf(x, y, Z);
    hold on
end
end


%% basis_function_ref
function N = basis_function_ref(x, y, basis_type, derivative_type)
if strcmp(basis_type, "P1")
    N = zeros(size(x,1), 3);
elseif strcmp(basis_type, "P1b")
    N = zeros(size(x,1), 4);
elseif strcmp(basis_type, "P2")
    N = zeros(size(x,1), 6);
elseif strcmp(basis_type, "P3")
    N = zeros(size(x,1), 10);
elseif strcmp(basis_type, "Non-com")
    N = zeros(size(x,1), 3);
elseif strcmp(basis_type, "Hermite")
    N = zeros(size(x,1), 10);
else
    error("Invalid basis type.");
end
if strcmp(basis_type, "P1")
    if strcmp(derivative_type, "x")
        N(:,1) = 1 - x - y;
        N(:,2) = x;
        N(:,3) = y;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1;
        N(:,2) = 1;
        N(:,3) = 0;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 1;
    end
elseif strcmp(basis_type, "P1b")
    if strcmp(derivative_type, "x")
        N(:,1) = 1 - x - y;
        N(:,2) = x;
        N(:,3) = y;
        N(:,4) = 27*x.*y.*N(:,1);
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -1;
        N(:,2) = 1;
        N(:,3) = 0;
        N(:,4) = -27*y.*(2*x + y - 1);
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -1;
        N(:,2) = 0;
        N(:,3) = 1;
        N(:,4) = -27*x.*(x + 2*y - 1);
    end
elseif strcmp(basis_type, "P2")
    lambda = 1 - x - y;
    if strcmp(derivative_type, "x")
        N(:,1) = lambda.*(2*lambda - 1);
        N(:,4) = 4*x.*lambda;
        N(:,2) = x.*(2*x - 1);
        N(:,5) = 4*x.*y;
        N(:,3) = y.*(2*y - 1);
        N(:,6) = 4*y.*lambda;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = 1 - 4*lambda;
        N(:,4) = 4*(lambda - x);
        N(:,2) = 4*x - 1;
        N(:,5) = 4*y;
        N(:,3) = 0;
        N(:,6) = -4*y;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = 1 - 4*lambda;
        N(:,4) = -4*x;
        N(:,2) = 0;
        N(:,5) = 4*x;
        N(:,3) = 4*y - 1;
        N(:,6) = 4*(lambda - y);
    end
elseif strcmp(basis_type, "P3")
    lambda = 1 - x - y;
    if strcmp(derivative_type, "x")
        N(:,1) = 0.5*lambda.*(-1 + 3*lambda).*(-2 + 3*lambda);
        N(:,4) = 4.5*lambda.*x.*(-1 + 3*lambda);
        N(:,5) = 4.5*lambda.*x.*(-1 + 3*x);
        N(:,2) = 0.5*x.*(-1 + 3*x).*(-2 + 3*x);
        N(:,6) = 4.5*x.*y.*(-1 + 3*x);
        N(:,7) = 4.5*x.*y.*(-1 + 3*y);
        N(:,3) = 0.5*y.*(-1 + 3*y).*(-2 + 3*y);
        N(:,8) = 4.5*lambda.*y.*(-1 + 3*y);
        N(:,9) = 4.5*lambda.*y.*(-1 + 3*lambda);
        N(:,10) = 27*x.*y.*lambda;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = -13.5*lambda.^2 + 9*lambda - 1;
        N(:,4) = 4.5*lambda.*(3*lambda - 6*x -1) + 4.5*x;
        N(:,5) = 4.5*x.*(6*lambda - 3*x + 1) - 4.5*lambda;
        N(:,2) = 13.5*x.^2 - 9*x + 1;
        N(:,6) = 4.5*y.*(6*x - 1);
        N(:,7) = 4.5*y.*(3*y - 1);
        N(:,3) = 0;
        N(:,8) = -4.5*y.*(3*y - 1);
        N(:,9) = -4.5*y.*(6*lambda - 1);
        N(:,10) = 27*y.*(lambda - x);
    elseif strcmp(derivative_type, "dy")
        N(:,1) = -13.5*lambda.^2 + 9*lambda - 1;
        N(:,4) = -4.5*x.*(6*lambda - 1);
        N(:,5) = -4.5*x.*(3*x - 1);
        N(:,2) = 0;
        N(:,6) = 4.5*x.*(3*x - 1);
        N(:,7) = 4.5*x.*(6*y - 1);
        N(:,3) = 13.5*y.^2 - 9*y + 1;
        N(:,8) = 4.5*y.*(6*lambda - 3*y + 1) - 4.5*lambda;
        N(:,9) = 4.5*lambda.*(3*lambda - 6*y - 1) + 4.5*y;
        N(:,10) = 27*x.*(lambda - y);
    end
elseif strcmp(basis_type, "Non-com")
    if strcmp(derivative_type, "x")
        N(:,1) = 2*x + 2*y - 1;
        N(:,2) = 1 - 2*x;
        N(:,3) = 1 - 2*y;
    elseif strcmp(derivative_type, "dx")
        N(:,1) = 2;
        N(:,2) = -2;
        N(:,3) = 0;
    elseif strcmp(derivative_type, "dy")
        N(:,1) = 2;
        N(:,2) = 0;
        N(:,3) = -2;
    end
elseif strcmp(basis_type, "Hermite")
    lambda = 1 - x - y;
    if strcmp(derivative_type, "x")
        a = x.*y.*lambda;
        N(:,1) = lambda.^2.*(3 - 2*lambda) - 7*a;
        N(:,2) = x.*lambda.^2 - a;
        N(:,3) = y.*lambda.^2 - a;
        N(:,4) = x.^2.*(3 - 2*x) - 7*a;
        N(:,5) = x.^2.*(x - 1) + 2*a;
        N(:,6) = x.^2.*y - a;
        N(:,7) = y.^2.*(3 - 2*y) - 7*a;
        N(:,8) = y.^2.*x - a;
        N(:,9) = y.^2.*(y - 1) + 2*a;
        N(:,10) = 27*a;
    elseif strcmp(derivative_type, "dx")
        b = y.*(lambda - x);
        N(:,1) = 6*lambda.*(lambda - 1) - 7*b;
        N(:,2) = lambda.*(lambda - 2*x) - b;
        N(:,3) = -2*lambda.*y - b;
        N(:,4) = 6*x.*(1 - x) - 7*b;
        N(:,5) = x.*(3*x - 2) + 2*b;
        N(:,6) = 2*x.*y - b;
        N(:,7) = -7*b;
        N(:,8) = y.^2 - b;
        N(:,9) = 2*b;
        N(:,10) = 27*b;
    elseif strcmp(derivative_type, "dy")
        c = x.*(lambda - y);
        N(:,1) = 6*lambda.*(lambda - 1) - 7*c;
        N(:,2) = -2*x.*lambda - c;
        N(:,3) = lambda.*(lambda - 2*y) - c;
        N(:,4) = -7*c;
        N(:,5) = 2*c;
        N(:,6) = x.^2 - c;
        N(:,7) = 6*y.*(1 - y) - 7*c;
        N(:,8) = 2*x.*y - c;
        N(:,9) = y.*(3*y - 2) + 2*c;
        N(:,10) = 27*c;
    end
end
end